function [ f, sensitivity ] = sensitivity_analysis( param, DATA, INFO, study_case, free, t_f1, t_f2, delta_t, var, LQL, activate_vd, use_Markov )

% INFO: function to do the sensitivity analysis

%% OUT:

    % f: (real) value of the perturbed cost function
    % sensitivity: (vec) sensitivity indices

%% IN:

    % param: (vec) initial parameters
    % DATA: experimental data (cell)
    %    DATA(:,1): days
    %    DATA(:,2): volume
    %    DATA(:,3): standard deviation
    % INFO: treatment info cell, each row represents a curve
    %    Columns: D; t_rad; anti-CTLA-4 IT (yes:1, no:0); t_treat_c4;
    %    anti-PD-(L)1 IT (yes:1, no:0); t_treat_p1
    % study_case: (vec) indices of the response curves that we want to study
    % free: (vec) controls free growth
    %   free(1): (int) 1 for free behavior until an specific time or volume, 0 otherwise 
    %   free(2): (int) 1 for time, 0 for volume
    %   free(3): (real) initial volume or time
    % t_f1: (real) time allowed to reach baseline tumor volume
    % t_f2: (real) simulation time
    % delta_t: (real) time step
    % var: (real) relative variation of the parameter
    % LQL: =1 if LQL, otherwise (modified) LQ (with PARAM(33)) 
    % activate_vd: =1 if vascular death (D > 15) 
    % use_Markov: =1 if Markov, otherwise deterministic

%% Algorithm

    f = zeros(length(param) + 1, 1);

    f(1) = least_squares(DATA, INFO, param, study_case, delta_t, free, t_f1, t_f2, LQL, activate_vd, use_Markov);
    h = param * var;
    param_var = param + h;
    sensitivity = zeros(1, length(param));
    param_in = param;
    for i = 1:length(param)
        param(i) = param_var(i);
        f(i+1) = least_squares(DATA, INFO, param, study_case, delta_t, free, t_f1, t_f2, LQL, activate_vd, use_Markov);
        sensitivity(i) = (f(i+1) - f(1)); %/(h(i));
        param = param_in;
    end

end


%%
function cost_tot = least_squares(DATA, INFO, param_new, study_case, delta_t, free, t_f1, t_f2, LQL, activate_vd, use_Markov)
    
    cost = zeros(1, length(study_case));
    par_c4 = param_new(25);
    par_p1 = param_new(35);

    for i = 1:length(study_case)
        j = study_case(i);
        D = cell2mat(INFO(j,1));
        t_rad = cell2mat(INFO(j,2));
        param_new(25) = cell2mat(INFO(j,3)) * par_c4;
        param_new(35) = cell2mat(INFO(j,5)) * par_p1;
        t_treat_c4 = cell2mat(INFO(j,4));
        t_treat_p1 = cell2mat(INFO(j,6));
        [sol, t_eq] = radioimmuno_response_model(param_new, delta_t, free, t_f1, t_f2, D, t_rad, t_treat_c4, t_treat_p1, LQL, activate_vd, use_Markov);

        if t_eq == -1
            cost_tot = 1e10;
%             disp('Initial tumor volume not achieved');
            return
        end

        ind = round(cell2mat(DATA(j,1)) ./ delta_t + t_eq / delta_t) + 1;
        data_y = cell2mat(DATA(j,2));
        err = cell2mat(DATA(j,3));
        sol = sol(ind);
        if ~isrow(data_y), data_y = data_y'; end
        if ~isrow(err), err = err'; end
        cost(i) = sum(((data_y - sol) .^ 2) ./ (err .^ 2));
    end
    
    cost_tot = sum(cost);
    
end